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Abstract. 

A computational model study for complete frequency redistribution linear 
incoherent two-level atomic radiation trapping in optically dense media using the 
multiple scattering representation is presented. This model study discuss at length 
the influence of the spectral distributions, overall opacity and emission quantum 
yield to trapping distorted ensemble quantities stressing physical insight and with a 
non-specialist audience in mind. Macroscopic reemission yield, lifetime, steady state 
spectra and spatial distributions arc calculated as a function of intrinsic emission yield, 
opacity and external excitation mode for Dopplcr, Lorentz and Voigt lineshapes. The 
work could constitute the basis for a final undergraduate or beginning graduate project 
in computational physics instruction and implements the analytical developments of 
the previous instalment of this contribution. 



PACS numbers: 02.50.Ey, 32.80.-t, 32.70.-n 



1. Introduction 



In optically thick media, electronic excitation energy can undergo several reabsorption 
and reemission events before either escaping to the exterior or being converted 
in thermal energy by means of collisional deactivation. This resonant radiation 
trapping is important in areas as diverse as stellar atmospheres (T], plasmas and 
atomic vapours luminescence [2], terrestrial atmosphere and ocean optics |3j, molecular 
luminescence [4j, infrared radiative transfer [5j and cold atoms [6]. 

The starting point of the majority of the incoherent radiation trapping models 
is the Holstein-Biberman equation which is a Boltzmann-type integro-differential 
equation describing the spatial and temporal evolution of the excited state number 
density. The previous instalment [7] of this contribution outlined the two alternative 
ansatze commonly used to obtain solutions for the classical trapping problem, 
Holstein's original exponential mode expansion and the so called Multiple Scattering 
Representation (MSR). The MSR solution was given a simple stochastic formulation 
and trapping dependent quantities (overall relaxation parameters such as ensemble 
emission yield and lifetime, time-resolved and steady-state spatial distributions as 
well as spectra) were calculated with the Holstein fundamental mode singled out. 
This instalment will now use a simple Markov chain algorithm to quantify incoherent 
trapping in a computational model study for two-level atomic models. It could 
easily be adapted into a computational physics project valuable in the context of 
atomic or computational physics instruction for final undergraduate and beginning 
graduate students. With this objective in mind, trapping dependent quantities are 
estimated in a unidimensional geometry for a single line with Doppler, Lorentz and 
Voigt spectral distributions. The excitation relaxation dynamics is considered for 
conditions mimicking electron impact as well as photoexcitation. A thorough discussion 
of the conditions for which the use of Holstein's fundamental mode alone is a tolerable 
approximation is included. 

Section [2] discusses the way the dynamics of incoherent trapping is taken into 
account within the framework of the MSR. It gives explicit expressions for the 
overall relaxation parameters as well as steady-state spectra and spatial distribution 
summarizing the previous instalment of this work. Section [3] explains the rationale as 
well as critical implementation details of the Markov stochastic simulation algorithm, 
the one chosen in this model study. Section [4] presents the results and their discussion at 
length. Particular emphasis is placed in the discussion of the physical implications of the 
spectral distribution as well as quantifying the relative contribution of the fundamental 
mode. Section [5] summarizes the main conclusions. Finally, the Appendices show 
implementation fine details and add some possible routes for complementing the work 
presented here. 
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2. Dynamics of incoherent trapping 



Radiation trapping studies should be cast in dimensionless coordinates since this 
increases computational efficiency and, more important, defines characteristic scales 
or universal conditions. The quantities most directly amenable to define characteristic 
scales in trapping are time, distance and optical frequency. The scaled time is t — Ft', 
where T is the global deactivation rate constant. The dimensionless distance is 
sometimes called the opacity or optical density and can be defined, along a given 
pathlength / and for homogeneously distributed species, as r = k (x) /$ (0) = 
k $ (x) /$(0). a; is used to represent the optical frequency (see below), k (x) is the 
single line monochromatic opacity (ko is the corresponding center-of-line value) and 
the absorption lineshape is given by the normalized spectral distribution <3> (x) (so 
that f*™&(x)dx = 1). For two-level atomic models the intrinsic, trapping 
undistorted, spectra can be written as a function of a dimensionless optical frequency 
defined as x = ^et- This is a normalized difference to the center of line frequency, 
where Au D stands for the FWHM of the Doppler distribution at each given temperature. 
For two-level atomic models, it is common to describe the absorption lineshape by 
Doppler - §d(x) = ~^ e ~ x2 ~ Lorentz - $l(^) = ~iq^2 ~ or Voigt - $y ( x ) = 

_|_ -u 2 

a 2 H !(. r „ M )2 d-u - spectral distributions. The Doppler distribution allow us to 
single out the pure Doppler broadening from the other broadening mechanisms while 
the Lorentz and Voigt 's distributions are the ones to be used in pure radiation damping 
and combined radiation and collision broadening conditions, respectively. In this last 
case, a = y^n (2)^^ is the Voigt characteristic width, the relative Lorentz over Doppler 
spectral width and implicitly dependent upon both temperature and vapour pressure 
via the dependence on collisional cross-section values. 

In the MSR ansatz for linear incoherent trapping with negligible time-of-flight for 
in-transit radiation, the spatial and temporal relaxation dynamics for excitation is given 
by 

n{r,t) = J2 a nPn(r)g n (t) , (1) 

n 

where n stands for the generation number of excited species (primordial excitation 
creates first generation, the trapping of this generation's reemission creates second 
generation and so forth; one can envisage each generation as the result of n — 1 
previous scattering events of resonant radiation), a n is the population efficiency for 
each generation, and p n (r) and g n (t) are the (normalized) spatial and temporal excited 
species distributions. 

From a practical point of view, the most important quantities are the macroscopic 
ensemble relaxation parameters, overall reemission yield and mean scaled lifetime r, 
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and the steady state emission spectrum I ss (x). These were derived in a previous 
contribution and are given by Pf): 



m ~ 1 a a 



gnfln + — , (2) 

n=l 1 a ™ c 

E^Ti 1 ng n a n + [m (1 - a nc ) + a nc ] 



and 



rSS (*) = z r VPi $ (*) • (4) 



J-T Sn=l # (*) ^ + 



$ fx) <ix 



In these, m stands for the first generation number that can be considered 
nonchanging (subscript nc is a remainder for nonchanging; see |7J) and a n = 
is the mean trapping or reabsorption probability. Finally, 



£ (*)=// e~^ r p n (r) drdS, (5) 
JqJv 

is the mean escape probability in the direction defined by the solid angle Q which, 
for a left escape from a unidimensional geometry (see next section), is just 

(fi(x)= [e-^ r p n (r)dr, (6) 



if one decides to start the opacity scale on the left side of the cell. 
It is also informative to know the steady-state spatial distribution, which can be 
cast as Ffl: 



nSS (r) _ En=l anVn (fj + J^Pnc (r) 



In the above expressions it is important to recognize that the trapping dynamics 
can be factored out into a generation varying part and another corresponding to 
generations that have the same spatial distribution (the fundamental mode) and the 
nonchanging part can be expressed explicitly as an analytical sum. The dynamics for 
this nonchanging part corresponds only to an attenuation of overall excitation in going 
from one generation to the next by a fixed a nc factor and, as a result, the contribution 
of this part corresponds to a monoexponential relaxation with a trapping dependent 
effective decay constant. 
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3. Markov stochastic algorithm 



The stochastic formulation of the Multiple Scattering Representation (MSR) allows to 
quantify the trapping influence on observables essentially trough the estimation of mean 
reabsorption and escape probabilities. In this respect, it is simpler and more amenable 
to discussion at an elementary level than the alternative Holstein multiexponential 
expansion. The MSR numerical implementation is straightforward and can be easily 
discussed at a level that stresses physical intuition without actually getting bogged down 
in all the technicalities of the troublesome estimation of Holstein modes other than the 
fundamental. We have furthermore decided to use an unidimensional geometry, driven 
by the motivation to do a simple mimic of a cylindrical tube of a discharge fluorescence 
lamp and by the desire to keep computational detail and power adequate to the proposed 
audience of final undergraduate and beginning graduate students. In this ID case the 
opacity scale comes naturally as the opacity along the cylindrical axis. 

The model used is a two-level single line atomic model with Doppler, Lorentz or 
Voigt spectral distributions and particular attention is paid to the influence of the 
lineshape on the trapping efficiency The Voigt case will illustrate the fact that a 
continuous variation of the characteristic width parameter will map the Doppler into 
the Lorentz distributions by changing the relative importance of the Lorentz-like wings 
over the Doppler-like core of the distribution (figure [I]). Accordingly, the ability to 
compute the Voigt line to machine precision is mandatory and Appendix A gives some 
implementation details in this respect. 

Complete frequency redistribution conditions are used [2], which means that 
the number of collisions during the lifetime of excited atoms is sufficiently high to 
render the reemitted photon's frequency completely uncorrelated with the frequency 
of the previously absorbed photon. In these conditions, the absorption and emission 
lineshapes coincide and the jump length distribution of the excitation random trajectory 
is independent of past history. This makes the formalism of Markov processes [8] 
especially adequate. Its rationale for the MSR implementation can be cast in the 
following way. The ID cell is divided into several bins, each corresponding to a pure 
state of the system and characterized by a mean probability that the excitation resides 
in that state. The system dynamics corresponds to the evolution of the probability 
of excitation being inside each state. The stochastic process is completely specified 
by (i) a column vector with the (normalized) spatial probability distribution of the first 
generation species, p\ = [p\], and (ii) a transition matrix P = [p 1 -?], whose entries are 
the one-step transition probabilities between states i and j. For complete frequency 
redistribution, there is an absence of memory effects (homogeneous chain) meaning that 
the transition probability between individual states depends only upon their relative 
opacity distance (it is independent of generation number and thus computed only once). 
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Figure 1. Doppler, Lorentz and Voigt lineshapes. 



The spatial distribution functions for all the generations are calculated from the previous 
generation by: 



The sample cell is divided into /i-length bins and the transition matrix elements 
are therefore given for a ID geometry by pi: 



corresponding to the Beer-Lambert law weighted according with the emission 
lineshape for an individual reemission-reabsorption (scattering) event. The integration 
takes into account all the possible emission frequencies, 1/2 is the left or right emission 
direction probability for a ID geometry and it was assumed that the bin width is 
controlled in order to attain a satisfactory precision. 

The complete specification of the Markov process is achieved once one specifies the 
initial spatial distribution p\. Two different cases were considered, one for homogeneous 
initial excitation (trivial) and another mimicking photoexcitation with the reabsorption 



p n+ i = Pp n . 





(9) 
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undistorted line. For photoexcitation from the left side of the ID cell, 

/+oo 
$ 2 (x)e-* {x)r >dx, (10) 
-oo 

which afterwards must be properly normalized. 

The binning of the spatial excitation distributions corresponds to the substitution 
of a continuous distribution for its discretized version, effectively transforming a 
continuous process into a discrete one realized in a lattice. It is therefore the Markov 



equivalent of a random walk defined over a regular spaced lattice 10 . The bin width (or 
the number of cells) is the critical parameter for the Markov algorithm. We have 
conducted several tests and found advisable to have a maximum bin size of 0.05 in an 
opacity scale. Otherwise, numerical artifacts associated with substituting the actual 
excitation migration for the jump between the mean coordinates of each bin of sample 
cell could exist. These were found to be the more important the higher the overall 
opacity. 

We have mentioned in the beginning of this section that the Markov 
implementation of the MSR model aims at estimating mean reabsorption and 
escape probabilities. The reabsorption probabilities are estimated with the following 
procedure. Each time ^ is used, the fraction of the excitation remaining inside sample 
cell gives the n th mean reabsorption probability a^. The excitation column vector is 
then (re)normalized and the process repeated. From the values of this parameter for 

ra-l 

each generation, the trapping population efficiency is = T7 a n- This procedure 

n=l 



is the equivalent of an importance sampling method (see Appendix B), in which it is 
assumed a unit intrinsic reemission yield in ([9]). The influence of the actual value of 
this yield (0o — f~+r — > ^ e ra ^° °f radiative over radiative plus nonradiative relaxation 
rate constants) is introduced analytically. Using the notation developed in the previous 
instalment, the Markov algorithm directly estimates a^, each generation population 
efficiency due to trapping (and geometry) alone, and the actual population efficiencies 
were then given by a n = a^0g _1 0- ^ s ^ or ^ ne esca P e probabilities, we have chosen the 
following implementation of (|6]). The monochromatic left escape probability is obtained 
from 

q^(x) = Q(x) Pn , (11) 

where p n is the spatial distribution and Q (x) is the escape matrix, whose entries 
are finally 

^(x)^e-*^, (12) 
which comes directly from Beer-Lambert law. 
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Equation (|8j) is equivalent to linear response theory and therefore the evolution of 
the spatial excitation is given as a convolution integral between the excitation profile 
p n with the delta response function given by ([£]). This renders the Markov approach 
more efficient since this convolution can be easily made using FFT algorithms [II] 
paying attention to zero pad to double size the column vector containing the excitation 



distribution in order to avoid wrap-around effects due to the cyclic convolution 12 



The speed-up factors could rise up to several orders of magnitude (roughly 50 to 400 
times for a number of Markov states of 2 000 to 200 000) making the FFT convolution 
the recommended implementation of the Markov algorithm. 



4. Results and discussion 



4-1. Ensemble relaxation 

Figure [2] shows the overall relaxation parameters for an initial homogeneous excitation 
and an intrinsic quantum yield of O = 0.9 as a function of opacity for Doppler, Voigt 
and Lorentz lineshapes. The higher the opacity the more important trapping is, with 
the following implications: (i) an increase of the mean relaxation time (equivalent to 
a mean number of scattering events before escape) and (ii) a decrease of the ensemble 
reemission yield (the fraction of original excitation that eventually comes out; an 
increased importance of trapping translates into additional possibilities of nonradiative 
relaxation). The Voigt continuous transition from Doppler into Lorentz is evident as 
well as the relative importance of core and wings of the distributions. The higher the 
Lorentz character of the spectra the higher the weight of the wings and the higher the 
escape probabilities (lower mean lifetime and higher macroscopic emission yield). 

Figures [3] and [4] show the mean scaled lifetime and reemission yield for homogeneous 
initial excitation for the Doppler and Lorentz limiting spectral distributions and for 
several values of the intrinsic quantum yield. Three main conclusions can be drawn 
from these results. First of all, the two most important parameters controlling trapping 
efficiency are the spectral distribution and the value of 0o- The higher the overall 
opacity the more difficult is the escape of radiation for Doppler-like distributions and 
the more important the escape from the Lorentz-like wings of the distribution. Second, 
trapping for Doppler-like distributions is much more efficient since, especially in high 
opacity cases, the escape of excitation at optical frequencies far enough from the line 
center frequency is reduced due to the extremely small probability of reemission at those 
frequencies (for unit reemission yield, the lifetime for the higher opacity is about 200 
for Doppler and only about 15 for Lorentz). Finally, under conditions rendering 
trapping efficient, the 0o value is of paramount importance; for unit intrinsic reemission 
probability all the excitation will eventually come out (thus giving a simple check for 
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Figure 2. Rccmission yield cj> and mean scaled lifetime r, with </>o = 0.90 for Doppler, 
Voigt a = 0.05, 0.1, 0.2 and 0.5 and Lorentz lineshapes (direction shown by arrow). 
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consistency of computation) but, as soon as 0o is smaller than one, each new scattering 
event gives the excitation another chance of thermal degradation (in ^ each generation 
contribution is intrinsically dependent upon 0o _1 )- Note in these figures that there is 
a strong dependence of the relaxation parameters with the 0o value, especially for the 
more trapping influenced Doppler case. 

All these conclusions are important in the discussion of atomic vapours ensembles 
for lighting applications, either electric discharge lamps or plasma display panel (PDP) 
devices. Better performance is achieved with higher macroscopic reemission yields. 
On top of that, the increase of the overall opacity is in principle desirable since this 
is related with the increase of the number of excited species. In a crude first order 
approximation, one can assume the overall lamp efficiency to be directly proportional 
to the product of (f) times the overall opacity (directly proportional to initial excitation 
density) : 

* oc 4> x r. (13) 

The actual behaviour of a lamp or a PDP can be quite involved since an increase 
in opacity means an increase of partial vapour pressure (the external dimensions of 
the device are fixed) and this induces several changes whose influence on the overall 
performance can be contradictory. The higher opacity means higher light throughput 
as long as the increase in the trapping efficiency does not substantially increase thermal 
degradation. The higher the opacity the more efficient the trapping (higher recapture 
probabilities) but, on the other hand, the resulting higher collisional effective rate 
constants could give rise to a smaller intrinsic 0o an d render the Lorentz distribution 
only approximately valid. A reduced </>o value and a shift of the Voigt distribution 
towards a more Lorentz-like profile tend by themselves to decrease trappping efficiency 
and increase light throughput. To have a simple idea of the effect, and due to the 
paramount importance of the reemission quantum yield (f>o, a series of results were made 
for both limiting Doppler and Lorentz distributions with the radiative quantum yield 
given by 0o = r r ^ T , where the quenching rate constant was assumed in a first order 
approximation to be linear with the cell opacity (T q = kr, with the numerical values 
T r = 10 7 s -1 and k = 2.5 x 10 4 s^ 1 ). This corresponds to the well known Stern- Volmer 
equation for dynamical quenching by binary collisions for unitary intrinsic radiative 
yield in the absence of collisions. The results are shown in figure [5] which is judged 
more representative of the actual lamp behaviour than the results in figures [3] and 
|4j Figure [5] shows that, ultimately, a delicate balance will dictate the best operation 
conditions which manifest themselves in the peaks of the \I/ values. Of course, the results 
are very approximate since the assumed functional dependence of 0o is only approximate 
and, even within the CFR two-level Lorentzian assumption, one should use an opacity 
dependent characteristic width of the Voigt distribution. Nevertheless, figure [5] allows 
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Figure 3. Rccmission yield <f> and mean scaled lifetime r for the Doppler lineshape, 
with 4>o values of 0.90, 0.93, 0.95, 0.98, and 1.0 (direction shown by arrow). 
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Figure 4. Rccmission yield (f> and moan scaled lifetime r for the Lorentz lineshape, 
with </>o values of 0.90, 0.93, 0.95, 0.98, and 1.0 (direction shown by arrow). 
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the discussion of the qualitative behaviour emphasizing physical insight without the 
additional burden of fine grained details. It shows how critical the spectral distribution 
shape and the quantum reemission yield are. For <p values sufficiently close to one, an 
increase in opacity corresponds to an increase in lighting efficiency due to the increased 
initial excitation number density. Due to the strong dependence of trapping on the O 
value, as soon as this value starts to be significantly smaller than one, trapping means 
an much increased importance of thermal degradation (compare the difference between 
the quantum and the ensemble yield or the reduced overall relaxation lifetime in the 
upper part of figure [5]). From some point onwards this will be more important than 
the increase in initial excitation due to a higher vapour pressure, which originates 
optimal operation conditions, giving the best possible lighting efficiency. Figure [5] also 
shows that the Doppler distribution is associated with smaller throughput in lighting 
applications when compared with the Lorentz case due to the step reduction of the 
ensemble reemission yield with the increase of the overall opacity. This is of course 
related with the use of an inert gas filling to render collisions more important (increasing 
the Lorentz character of the spectral distribution and reducing trapping efficiency) in 
fluorescence lamps. 

Finally, figure [6] shows the predicted ensemble relaxation parameters for both 
homogeneous and photoexcitation as a function of overall opacity. Up to opacities of 
the order of 10 no significant difference exists (photoexcitation is able to penetrate well 
deep into sample cell). But, for higher opacities, the importance of trapping continues 
to increase indefinitely for homogeneus excitation while it levels off for photoexcitation, 
a point to be revisited in section |4.3| when discussing spatial distribution functions. 



4-2. Steady-state spectra 

Figure [7] shows the estimated normalized spectral distribution in steady-state conditions 
for Doppler, Lorentz and Voigt lineshapes for both primary homogeneous and 
photoexcitation. The motivation for the homogeneous case is the excitation along 
the axis of a fluorescence lamp for lighting applications. It shows the well known self- 
reversal of spectral lines due to the higher attenuation of core optical frequencies. For 
photoexcitation there is a balance between reduced penetration of external excitation 
and higher attenuation at core frequencies which dictates a flattening of the spectra near 
the line center (of course, for left wall photoexcitation and right wall detection there is a 
self-reversal higher than the one for homogeneous excitation; not shown). In both cases, 
there is a considerable broadening of the detected spectra and the Voigt distribution 
has an intermediate character between core Doppler-like and wings Lorentz-like. 
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Figure 5. Ensemble relaxation parameters as a function of overall opacity for Doppler 
and Lorentz lineshapes. The upper part shows the Reemission yield <p and mean scaled 
lifetime r while the lower part shows the predicted relative efficiency for lighting 



applications estimated from (13). </>o is implicitly dependent upon opacity (upper 
part and text). 



14 



1 I 1 — I — I I I 1 1 II 1 — I — I I I 1 1 II 1 — I — I I 1 1 1 II 1 — I — I I 1 1 1 1 

AAAAA, ' ' ' " 



0.1 - 



A A, 



Doppler 



a a Aaaa aaaaaaaa aaaaaaaaaaa 



A 



A 



Photoexcitation 



A 



Homogeneous 
excitation 



A 



A 



A, 



A, 



A. 



' A AAA 



-i i i Mini i i i Mini i i i Mini i i i i i 1 1 1 



10° 



l(f 



10 
8 - 

6 - 
4 - 



2 - 



10 1 

~ I — I | 1 1 — I | 1 | 1 l . \ M/ 



10 4 



xAA 



aA 



aa 



AAA 



aA Aaaa™ 



Homogeneous 
excitation 



A 



A 



A. 



a a aAaaaaaaaaaaaaa aaaaaa aa aa| 



Photoexcitation 



A 



A 



-J i i Mini i i i Mini i i i Mini i i i ii 1 1 1 



10° 



10 1 



Opacity 



10 d 



10 4 



Figure 6. Rccmission yield cf> and mean scaled lifetime r for the Doppler lincshapc 
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Figure 7. Normalized steady-state spectra for Doppler, Lorentz, and a = 0.1 
Voigt lineshapes for an overall opacity of 100 and <po = 0.90. The photoexcitation 
case corresponds to both excitation and detection from the left cell wall using for 
photoexcitation the reabsorption undistorted line. 
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4-3. Spatial distribution 

Figure [8] shows the spatial distributions for both limiting cases of Doppler and 
Lorentz distributions. This figure shows that the fundamental mode spatial 
distribution (limiting case for a relaxed, non-changing spatial distribution and thus 
independent upon the original excitation) could only give a reasonable approximation 
of the steady-state distribution for the homogeneous excitation case; for photoexcitation 
it is more convenient to choose the spatial distribution of the first generation species 
as a first approximation to the overall distribution. This illustrates the well known 
procedure of, whenever approximating the actual trapping dependent behaviour by 
the monoexponenial fundamental mode (easier to obtain by a variational procedure 
or given by Holstein's asymptotic approximations), design the experimental setup to 
mimic as much as possible the fundamental mode spatial distribution (symmetrical 
and well spread into the bulk of sample cell) with the external excitation. This can 
be accomplished with photoexcitation of high opacity samples using strongly detuned 
external radiation. 

Figure [8] shows some of the difficulties of quantifying trapping simply by using the 
fundamental mode, a point further illustrated in table [I] Several conclusions can be 
draw from its data: (i) Doppler distributions render trapping much more efficient and 
thus its fundamental mode contribution is always much higher than for the Lorentz case, 
(ii) the spatial spreading for Doppler is smaller, giving rise to higher generation number 
for the fundamental mode, (iii) the use of the fundamental mode alone for Lorentz 
distributions (and therefore, albeit with a lesser degree, for Voigt) is never justified and 
(iv) to approximate the actual behaviour for photoexcitation to the fundamental mode 
is never justifiable. 

Two additional points related with the common practice of using only the 
fundamental mode to take into account trapping distortions should be stressed out. 
First, the fundamental mode is the slowest decaying possible and is located well (and 
symmetrical) into sample cell. To substitute the whole of the ensemble dynamics for 
the fundamental mode alone will always overestimate the lifetime and underestimate 
the reemission yield (spatial distribution giving the highest possible trapping efficiency) 
thus introducing a systematic error. Secondly, the use of the fundamental mode alone 
is too often misunderstood with the use of the asymptotic approximations proposed 



by Holstein 13 , only valid in the high opacity limit and for ideal geometries [2j. The 
MSR has a clear cut advantage in this respect since it allows an easy estimation of 
the fundamental mode as the one corresponding to a nonchanging spatial distribution, 
irrespective of opacity and geometry. 

The spatial distribution functions presented in this section draw some further 
insight into the previous results of figure [6j The leveling of the lifetime and reemission 
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Figure 8. Normalized spatial distributions of excitation in steady-state (SS) 
conditions for Dopplcr and Lorentz lincshapes for an overall opacity of 250 and 4>q = 
0.90. It is also shown the primary excitation (homogeneous or photoexcitation) as 
well as the fundamental mode distributions. The photoexcitation case corresponds 
to both excitation and detection from the left cell wall using for photoexcitation the 
reabsorption undistorted line. 
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Table 1. Fundamental mode contribution to the reemission yield <j> and mean scaled 
lifetime r. Also shown is the approximate generation number corresponding to the 
fundamental mode (m), for a 10 -6 fractional tolerance to consider a non-changing 
spatial distribution (see text). In all cases, <fi = 1. 

Homogeneous Photoexcitation 

Opacity Doppler Lorentz Doppler Lorentz 

m(pTm(f)Tm(f) r m r 

10 10 6% 20% 7 1% 3% 10 3% 14% 7 1% 3% 
100 70 10% 40% 30 0.5% 3% 100 2% 20% 35 0.1% 1% 
1000 200 25% 60% 80 1% 6% 550 0.2% 6% 100 0.1% 1% 



yield at higher opacities for the photoabsorption case corresponds to absorption of 
external excitation complete within a layer smaller than the overall opacity making 
the ensemble relaxation effectively insensitive to the overall opacity. This approaches 
well the conditions of semi-infinite geometry, under the time scale of ensemble complete 
relaxation. 

5. Conclusions 

This instalment presented a unidimensional computational model study for complete 
frequency redistribution linear incoherent atomic radiation trapping illustrating the 
numerical implementation of the stochastic model developed previously. It illustrates 
the advantages of the multiple scattering representation (MSR) over the Holstein 
expansion, based on physical insight and computation feasibility at an elementary 
level. Holstein's ansatz has significant shortcomings when compared to the equivalent 
alternative of MSR: (i) Holstein spatial modes are unphysical except for the 
fundamental, (ii) their estimation is computationally much more troublesome than 
the simple algorithms used in this work, (iii) the wide spread use of original Holstein 
expressions for the fundamental mode are only valid in the asymptotic limit of high 
opacities while MSR allow an easy estimation of the fundamental at any opacity value 
and (iv) the higher Holstein modes are difficult to obtain while for Lorentz-like spectral 
distributions we found that their contribution must be always taken into account (the 
fundamental mode contribution to ensemble relaxation being always small; higher 
Holstein modes correspond in the MSR language to small number generations and 
are easy to obtain with the stochastic formulation presented). 

The dependence of the ensemble reemission yield and lifetime, relative efficiency 
for lighting applications, steady-state spectral and spatial distributions on quantum 
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yield, opacity and homogeneous or external photoexcitation are discussed at length 
for the Doppler, Voigt and Lorentz lineshapes. We quantify the contribution of the 
nonchanging fundamental mode and found troublesome using uniquely this mode for 
Voigt and Lorentz like spectra. The results should appeal to a broad audience and 
provide insight in a wide range of more realistic situations in an atomic as well as in an 
astrophysical context. 

Several possible developments of the general framework presented in this work are 
shown in Appendices B to D. In the first, we give an outlook of the implementation 
details that we have found critical for a Monte Carlo simulation alternative of the 
Markov algorithm. The Monte Carlo and Markov approaches constitute in fact two 
general purpose algorithms for incoherent radiation propagation problems. Each has in 
own advantages and shortcomings. The Monte Carlo constitute a simulation of particle 
like trajectories while the Markov model quantifies the evolution of mean probabilities. 
We advocate the second alternative in all but the more demanding cases (detailed 
3D geometries, realistic multi-level atomic models, partial frequency redistribution and 



polarization dependent radiation transport, issues not addressed here). In Appendix 



[Cj we give a brief description of a more realistic but nevertheless still unidimensional 
Markov implementation of a radiation transport model, especially appealing for plane- 
parallel stellar atmosphere theories. Finally, in Appendix D we give an analytical 
modification of the Markov algorithm for multiple specular reflections in boundaries. 
The appendices could be useful for more advanced projects. 
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Appendix A. Numerical Voigt distribution 

The numerical evaluation of the Voigt spectral distribution can be troublesome, as one 
can easily judge from the large number of approximations that have been published in 
the last decades balancing precision and computational speed (see [2] and references 
therein). This is especially true for the wings of the distribution in case of trapping 
since the photons can most easily escape via the wings, especially in high opacity 
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vapours. However, given the current desktop computing capabilities, the numerical 
(careful) direct integration of the defining equation is perfectly adequate and the use of 
approximations to reduce computation time is not justifiable any more. The difficulty in 
the direct numerical integration cames from the behaviour of the integrand function: it 
differs from zero over two width scales, a broad scale centered in zero (corresponding to 
the exponential term in the numerator of the integrand) and a much narrower one 
centered in a frequency corresponding to the Voigt frequency (associated with the 
difference term in the denominator). The integration domain should thus be broken 
into smaller domains and an automatic adaptative integration algorithm should be 
used in each subdomain always starting at the integrand function maximum and with 



a initial stepsize adapted to the local scale of variation of the integrand 12 . We 
have used the 400 central frequencies for the central broad feature and a 0.4 frequency 
width for the floating narrow peak. Integrations further away from the central core 
were analytically mapped from an infinite to a finite integration range. In order to 
decrease the time for repetitive Voigt functions evaluations, the Voigt distribution was 
previously computed in a given table of frequency values and, whenever necessary, cubic 
spline interpolated [12]. We used a linear scale in the core (frequency range up to 100 
with a 5 x 10 -2 spacing) and a log scale in the wings (frequency range from 100 to 
10 8 with 1 x 10~ 3 logxo spacing). Natural cubic spline was not necessary since the 
derivatives of the Voigt distribution in the end points are analytical. 

Appendix B. Monte Carlo simulation 

The mean reabsorption and escape probabilities can be estimated either with a Markov 
chain algorithm or with Monte Carlo simulation mimic of the experiment. We will 
outline the computational details we found critical for the Monte Carlo alternative and 
give our best advice on each method strengths and limitations and under which physical 
conditions the Monte Carlo is especially adequate. 

The Monte Carlo (MC) method makes a direct simulation of the trapping process 



using particle-like trajectories for radiation in cell 14 . The initial excitation coordinate 
is randomly chosen corresponding to either homogeneous or external photoexcitation. 
The (re)emission coordinate is the same as the absorption one and, after emission, 
a random direction and optical frequency x must be chosen from the appropriate 
distributions. The photon path in cell is followed and a reabsorption coordinate 
drawn from the Beer-Lambert exponential distribution with absorption coefficient given 
by <£>(x). This should be then tested for escape from cell; if the photon escaped, 
another excitation trajectory should be initiated from the first generation, otherwise the 
simulation must continue by increasing generation number and repeating. Appropriate 
counters keep track of the actual number of trajectories giving rise to n th generation 
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species (the ratio to the total number of trajectories initiated for first generation gives 
the population efficiencies), and the mean escape probabilities (for each generation 
number, excitation coordinate r and optical frequency x, the escape probabilities 
counters are incremented with e _ *^ r for left escape; compare (JfjJ). To decrease the 
computation burden, an importance sampling method should be used always assuming 
a unit intrinsic reemission yield and after the influence of the actual value of the yield 
introduced analytically, as described in the main text. The spatial distribution functions 
for each generation are most easily obtained for a discretized cell just by keeping track 
of the number of excited species in each cell bin but it is important to acknowledge 
the fact that this binning is used only for graphical representation purposes since the 
mean escape probabilities are computed from the actual spatial coordinates. This is an 
an important difference relative to the Markov algorithm since this last case effectively 
makes a simulation on a lattice model from the very onset. 

The complete Monte Carlo description lacks only the fine details of the 



implementation of the transformation method to obtain non-uniform deviates 12 



These are used in the simulation to obtain excitation coordinates and optical 
frequencies. The excitation coordinates are either homogeneous for initial 
excitation (trivial) or drawn form the Beer-Lambert law, giving the reabsorption 
coordinate relative to emission point in an opacity scale for a given absorption coefficient 
at each optical frequency Q(x). Since Beer law corresponds to an exponential 
distribution, the random deviates are analytically given by —ln(y)/Q (x), where y 
is a uniform deviate. For the optical frequencies x, we must distinguish between 
Doppler, Lorentz and the Voigt lineshapes. The Doppler and Lorentz distributions 
have analytical inverses and are therefore simpler. Lorentz deviates are given by x = 
tan \k (y — 1/2)] while Doppler correspond to zero mean 1/2 variance normal deviates. 
The reference [12] only gives directly a routine for zero mean unit variance deviates 
but it is easy to derive a general purpose formula for zero mean, arbitrary a variance 
along the same lines of the used Box-Mueller algorithm. The Voigt case is more 
troublesome since the transformation method must be implemented numerically. The 
most straightforward implementation of this is to use a cubic spline of the cumulative 
distribution function, truncated to a very high frequency (se have obtained good results 
with an upper limit of 10 8 , on a logarithmic scale in the wings). 

For the simple model study of this work, the Markov algorithm has distinctive 
advantages over the Monte Carlo alternative since it considers directly the evolution 
of mean excitation probabilities. These can be obtained with a fast algorithm and 
this is very important for trapping due to the critical need to proper quantify the 
fundamental mode. The estimation of the fundamental mode in MC is much more 
difficult due to the characteristic slow convergence of MC estimates and to the fact 
that MC simulations must be made with a maximum generation number specified at 
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start. This is especially important for high opacity cases since in this case the MC 
computation time can be prohibitly large and we have found that a small error in 
the fundamental mode can be greatly amplified in overall relaxation parameters. MC 
is nevertheless more versatile with respect to generalizations to either more realistic 
geometries and atomic models as well as to include polarization effects. The influence 
of geometry deserves a small note. For realistic tridimensional geometries the demands 
introduced by a sufficiently fine spatial discretization in the Markov procedure can 
compromise the practical feasibility with current desktop technology. For MC there is 
no significant added overhead computational complexity since the level of discretization 
only changes the visual resolution of the spatial distribution functions. Balancing the 
shortcomings of both approaches, our advice would be to use Markov whenever possible. 
For 3D geometries, a possible procedure would be: (i) to identify the smallest dimension 
in cell, (ii) to do a ID Markov estimation of the fundamental mode generation number 
for that smallest opacity and (iii) finalize by MC simulating the actual geometry until 
that generation number (assumed to describe well the fundamental mode). 



Appendix C. Unidimensional models for radiation transport 

The unidimensional implementation of radiation migration used in this work is certainly 
very naive. An alternative ID possibility nevertheless exists that does not imply a 
qualitative increase in computation complexity and which is more realistic. It introduces 
some approximations, albeit of a different nature. To keep the unidimensional 
formulation of radiation migration, one can represent trapping as a function of a 
characteristic distance measured in the perpendicular to some known surface. This 
corresponds to the well known cases of plane parallel stratified stellar atmospheres [l] 
or the idealized one-dimensional geometries in laboratory scale atomic vapours pi. In 
the Markov algorithm, the one-step transition probabilities should be modified in order 
to take into account all the possibilities of transition between any two coordinates with 
the same difference in perpendicular distances to the reference surface. This can be done 
for homogeneous three-dimensional space, by projecting the transition probability on 
an arbitrary axis (the one used for the ID geometry). Instead of ([9]), one obtains 15 

1 r+oo 

P v ~ -h / $ 2 (x) E 1 (|r< - Tj\ $ (x)) dx, (C.l) 

where Ex (x) is the exponential integral function, defined as Ei (x) = f^ 00 — du. 
No further modifications of the Markov stochastic algorithm are needed. 
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Appendix D. Specularly reflecting walls 



Another case of practical interest is a vapour cell with partially reflecting walls (2j. For 
the unidimensional geometry and for specularly reflecting walls, the Markov algorithm 
is amenable to a simple modification since the absorption, transition and escape 



probabilities can be cast as series which have analytical representations 16 



Consider a unidimensional geometry with zero at the left and r max maximum overall 
opacity in the right wall. The corresponding reflectances are assumed constant and will 
be represented as Rl and Rr. The stochastic interpretation of the Beer-Lambert law 
gives for the absorption and escape mean probabilities, after an r optical pathlength, 
$ (x) e~®( x ' r and e~*( x ) r , respectively. $ (x) is just the mean probability of reemission 
of a photon with frequency x. With all of this in mind, the Markov transition and escape 
probabilities can be reformulated taking into account the possibility of multiple specular 
reflections of the geometry boundaries. For the mean transition between Markov states i 
and j, instead of (|9|) one can write 



I r+oo 

P v~-h r>(x,R L ,R R )dx, (D.l) 
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with the integrand function given by 
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The first line gives the transition probability due to direct absorption prior to 
any reflection. The other lines give the additional absorption probability in state j 
after at least one reflection in the cell walls. The second and third lines correspond to 
emission from state i to the left while the fourth and fifth correspond to emission to the 
right. Finally, the second and fifth lines consider photon absorption in state j with the 
radiation coming from the left of this state while the third and fourth lines correspond 
to absorption from the right. 

Rearranging, one eventually obtains 
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P {x,R L ,R R ) = & ( x )e-"^-^ + 
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Using the same procedure, the left escape probability is given by 
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The first line is the Ze/t escape, given that the emission was is state i to the left 
while the second is the left escape for initial right emission. 

Finally, multiple reflections can also change the initial spatial distribution in the 
case of external photoexcitation. For the case of photoexcitation from the left side with 



the undistorted resonance line, (10) should be modified into 
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The first line is the absorption in state i of radiation propagating from left while 
the second is the absorption of radiation coming from the right. 

For specular reflection in this unidimensional model the advantages of the Markov 
algorithm over the alternative Monte Carlo are even more important than for the case 
in which reflection in boundaries is not considered. The multiple reflections are taken 
into account analytically for the Markov case while in the Monte Carlo simulation the 
increase in computation time with the increasing importance of reflections can be very 
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high. Nevertheless, the increase in trapping efficiency with the added possibility of 
specular reflection is only significant if both reflection coefficients are important. In 
fact, a ID geometry in which one of the boundaries is perfectly reflective while the 
other has a zero reflection coefficient is equivalent to a ID sample of twice the size of 
the original one. 
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